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Abstract 

We study the 0(N) loop model on the Honeycomb lattice with real value N > 1 
by means of a cluster algorithm. The formulation of the algorithm is based on the 
equivalence of the 0(N) loop model and the low-temperature graphical representa- 
tion of a iV-color Ashkin- Teller model on the triangular lattice. The latter model 
with integer N can be simulated by means of an embedding Swendsen- Wang- type 
cluster method. By taking into account the symmetry among loops of different col- 
ors, we develop another version of the Swendsen- Wang- type method. This version 
allows the number of colors N to take any real value N > 1. As an application, 



* Correspondence should be sent to: ydl0@nyu.edu 



we investigate the N = 1.25, 1.50, 1.75, and 2 loop model at criticality. The deter- 
mined values of various critical exponents are in excellent agreement with theoretical 
predictions. In particular, from quantities associated with half of the loops, we de- 
termine some critical exponents that corresponds to those for the tricritical q = N 2 
Potts model but have not been observed yet. Dynamic scaling behavior of the 
algorithm is also analyzed. The numerical data strongly suggest that our cluster 
algorithm hardly suffers from critical slowing down. 

1 Introduction 

In Monte Carlo studies of statistical systems undergoing phase transitions, critical 
slowing down is one of the prominent problems. Consider a Monte Carlo algorithm with 
dynamic exponent z > 0. In order to generate a given number of effectively independent 
samples, one has to spend computing effort oc L d+Z , where L is the linear sytem size and L d 
accounts for the volume of the sytem of interest. For the local Metropolis algorithm for the 
Potts model, the dynamic exponent is around 2.2. Thus, in two dimensions, the required 
computing effort grows like L 4,2 . Therefore, a central task of computational statistical 
physics is to develop algorithms such that z vanishes or is significantly suppressed. For a 
discussion of Monte Carlo methods, see Ref. [1]. 

For the Potts model, a significant breakthrough was the invention of the Swendsen- 
Wang cluster method [2] and its single-cluster version-the Wolff cluster method [3]. For 
the two- and three-dimensional Ising model, the dynamic exponent of the SW algorithm 
is about 0.1 and 0.45 [4], respectively. In comparison with the Metropolis simulations, 
the critical slowing down is significantly suppressed. 

In addition to the Potts model, another important class of models in statistical physics 
is the O(N) model. The O(N) model is denned in terms of iV-component spins on 
a lattice, with an isotropic pair coupling of the form = e(Si ■ Sj), where % and j 
are a pair of neighboring lattice sites and e is a function. A particularly interesting 
case is the honeycomb O(N) model, where function e is e(p) = — log(l + xp), with x a 
measure of the inverse temperature. It turns out that the O(N) model has a nice graph 
representation [5-8]; the graph consists of a number of nointeracting and nonoverlapping 
loops on the honeycomb lattice. 

However, in contrast to the Potts model, an efficient cluster algorithm is still lacking 
for the O(N) model. When simulating the O(N) model, one has to apply a Metropolis- 
like local algorithm, except for some special cases such as N = or iV = 1. Even worse 
is that local updates of loop configurations require some global connectivity information. 
Therefore, the computing effort grows like L d+Z+Z ' as L increases, where the exponent z' 
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accounts for the effective critical slowing down due to the global-connectivity-checking 
procedure. The value of z' is close to 2 in two dimensions, unless some complicated data 
structure is applied. 

Apart from the above computational considerations, developing efficient algorithms for 
the O(N) loop model is also highly desirable from physical point of view. This is because, 
while much exact information about the critical properties of the O(N) loop model has 
been accumulated in the past decades, many open questions still exist; for some recent 
publications, see e.g., Refs. [9-13]. As a generalization, dilution can be introduced into 
the O(N) loop model and the Potts model. For a sufficient number of diluted sites, a 
different universality class can arise at the so-called tricritical point. In contrast to the 
tricritical Potts model, exact results for the tricritical O(N) model are scarce. Therefore, 
a high-precision numerical study of the O(N) model is of great importance. 

In this work, we solved the long-standing problem of developing an efficient algorithm 
for the O(N) loop model for N > 1. We apply the newly developed cluster algorithm, 
an embedding Swendsen- Wang- type cluster method, to the honeycomb O(N) loop model 
without dilution, for which many exact predictions are available. The numerical data 
confirm the exact predictions by Coulomb gas theory and by conformal field theory. Fur- 
ther, the dynamical data imply that, for N > 1, the embedding cluster algorithm hardly 
suffers from critical slowing down. 

2 The O(N) loop and the Ashkin- Teller model: exact 
mapping 

Let us consider a plane graph Q = (V, E) and its dual graph Q* = (V*,E*). Let |V| 
and \E\ denote the total numbers of vertices and of edges of Q, respectively, and \V*\ and 
\E*\ for graph Q*. The dual relation between Q and Q* guarantees \E\ = \E*\; there is 
one-to-one correspondence between edges of Q and Q* . 

On the edges of Q, bonds are placed such that they form a number of closed paths, 
or loops. In addition, it is required that these loops cannot share a common bond 1 . The 
weight of such a loop configuration is given by w l N c , and the partition sum is 



where the sum is over all satisfied loop configurations. Symbol £ denotes the sum of the 



Loops are allowed to intersect in the sense that they can share a common vertex. 




(2.1) 



loops 
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lengths of all loops, and c is the number of loops. In principle, parameters w and N can 
be any real or complex numbers. In this work, we shall only consider real numbers w > 
and N > 0, so that we have a probabilistic interpretation. 

First, let us further restrict our attention to the case where N > 1 is an integer, and 
define a iV-color Ashkin- Teller (AT) model [14, 15] on the dual graph Q* 2 . On every 
vertex i E V* of Q*, one simultaneously places N independent Ising spins say they 
have color m — 1, 2, • • • , N. On every edge (i, ) 6 E*, spins in the same color interact via 
couplings J2, and any two pairs of spins and that are in different colors interact 
via couplings J4. The latter involves four-spin interactions. The Hamiltonian reads 

N 

n AT (J2,j„N) = -j 2 j2 E *< (m M m) - J *E E *l m M m M B M B) • ( 2 - 2 ) 

m=l {i,)eE* m>n (i.)e_B* 

The partition sum is then given by 

^at(j 2 , j 4 ,ao = E e_WAT(<T) ' ( 2 - 3 ) 

where the sum is over all spin configurations cr*™-* with m — 1, 2, - • • , N, as represented 
by a single symbol cr. Apparently, in the low-temperature region J 2 , J 4 — > 00, systems 
defined by Eq. (|2.2j) have A^ 2 ground states. The N = 2 case reduces to the isotropic 
version of the standard Ashkin- Teller model [14,15]. 

In this work, we shall concentrate on model ()2.2|) in the infinite-coupling limit: J 2 —* 
00, J 4 — > —00, but J2 + (N — 1) J4 = J fixed at a finite value. In order to see the physical 
implication of this limit, let us consider the spin configurations on the ends of a given 
edge e of graph Q* . From Eq. (|2.2|) . if there are k unequal Ising variables in the same 
color on edge e, the Boltzmann weight reads 

k = : exp[J 2 N +^J 4 N{N -1)} 

k = l : exp{J 2 (AT-2) + ij 4 [A r (A r -l)-2(Ar-l)]} 

fc : exp{J 2 (A^-2A;) + ^J4[A r (A r -l)-2A;(A^-l) + A;(A;-l)]} . (2.4) 

Normalizing these weights by that of case k = 0, one has Tabled 

2 In Refs. [14,15], only N = 2 was considered. Nevertheless, the Ashkin- Teller model with N 7^ 2 has 
also received considerable research attentions; see e.g., Refs. [16,17] 
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Configuration 


Normalized weight 


k = 


1 


k = 1 


e -2[J 2 + (A r -l)J 4 ] = e -2J 


k > 2 


g— 2kJ _ gfc(fc— 1)J4 y q 







Table 1: Boltzmann weights for the model (|2.2j) in the infinite-coupling limit. For any 
chosen two pairs of spins and cr^ n ^ with m^nwe first quote the weight read off from 
([2.2)1 . and then the normalized weight obtained by making the first one equal to 1. We 
have used J to specify J 2 + (N — 1) J4. Notice that the weight for k > 2 vanishes. 

In words, for any given edge e G •£"*, only two types of spin configurations are allowed 
on its ends: spins in the same color and are all equal, or there is at most one 
pair of unequal spins in the same color. The relative weight of the latter over the former 
is e~ 2J . From now on, we shall refer to the AT model whose relative weights are given in 
Table Qas the infinite-coupling AT model (IAT). 

It turns out that, at least for integer N, the loop model defined by Eq. (|2.1|) can be 
exactly mapped onto the iV-color IAT model. For this purpose, let us consider the low- 
temperature expansion of the IAT model, which can be performed in a similar way to the 
low-temperature expansion of Onsager's Ising model (see e.g., ref. [18]). For each color of 
Ising spin variable a( m \ we represent those unequal neighboring spins by lines of color m 
on the corresponding edge of the graph Q (recall that the IAT model itself is defined on 
the dual graph Q*). In words, if two adjacent spins of color m are unequal, draw a line of 
color m on the edge of Q; otherwise, do nothing. Do this for all pairs of nearest-neighbor 
spins. Apart from effects caused by boundaries 3 , one obtains an Eulerian subgraph E'\ 
i.e., for every vertex, the number of lines touching it must be even. The lines are therefore 
joined together to form polygons (loops); these polygons have color m. 

Conversely, these polygons divide the plane into spin-up and spin-down domains for 
Ising variable a^ m \ For any such set of loops, there are just two corresponding spin 
configurations: they are related to each other by flipping all spins a^ m \ 

Do this for all colors of Ising variables a (recall that we have used symbol a to represent 
all spin variables for m = 1, 2, • • • , N). One obtains a configuration containing TV- 
color loops. Any such set of loops corresponds to iV 2 spin configurations. 

For the IAT model, the zero weight of configuration 7^ aj"^ , 7^ crj™^ for any 
m 7^ n guarantees that there is at most one pair of unequal spins in same color on each 

3 One would expect that boundary effects vanish in the thermodynamic limit. 
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edge of Q*. Therefore, although loops on graph Q can intersect, they can never share a 
common edge. An example of the low-temperature graph for the N = 2 IAT model on 
the triangular lattice is shown in Fig. ^ 




: (T (D = _ lj(7 (2) = _ 1 | . a (i) = ha (2) = _ 1 | (J (i) = _ 1)(7 (2) = 1 



Figure 1: A low-temperature graph of the N = 2 IAT model on a triangular lattice of 
size 6x6. The edges of the dual lattice - the honeycomb lattice, are shown as solid gray 
lines. Loops corresponding to Ising variables and are displayed as thick solid 
lines in red and in blue, respectively. No loops can intersect on the honeycomb lattice. 



Let the energy of a ground state be zero. Table tells that, in comparison with a 
ground state, a pair of unequal spins must receive an energy penalty e~ 2J , denoted as w. 
The partition sum of the IAT model can then be written as 

loops {-r,} 

where the second sum is over all possible color arrangements for a given loop configuration. 
The factor N 2 accounts for the number of ground states. Note that each loop can randomly 
take one of the N colors with equal probability, and the second sum can then be easily 
evaluated. This leads to 

Z AT (J,N) = iV 2 ^w/iV c , (2.6) 

loops 
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where £ and c again represent the sum of the lengths of loops and the number of loops, 
respectively. 

Apart from a trivial constant, the partition sum (|2.fi|) is equal to Eq. 1)2.1)1 . In other 
words, an exact one-to-iV 2 mapping has been established between the loop model (j2.1j) 
and the iV-color IAT model defined by Table H The weight w of the line segment in the 
loop configuration is related to J by w = e~ 2J . 

Note that, if the colors of loops are ignored, the loop representation of the IAT model 
is also a low-temperature graph of the Ising variable s = Yi m cr< ' m ' ) > ^he P r °duct of all N 
colors of variables . This observation is vital in the development of cluster algorithms 
for the O(N) loop model, as shown later. 

3 Cluster simulation of the AT model 

Given the well-known Swendsen-Wang or Wolff cluster algorithm for Onsager's Ising 
model, an analogous cluster algorithm is readily available for the AT model described 
by Eq. ()2.2)1 4 . The so-called direct or embedding algorithms can be formulated. A 
detailed stuty of the Swendsen- Wang-type algorithms, including its dynamical behavior, 
was carried out by Salas and Sokal [20] in the context of the N = 2 AT model on the 
square lattice (not in the infinite-coupling limit). In the present work, we shall only 
consider the embedding version of the Swendsen- Wang-type algorithm. 

For the AT model ()2.2jl . for a given color of Ising variable a^ m \ the effective Hamilto- 
nian, conditioned on the other spin configurations for n — 1, 2, m — 1, m + 1, ■ • • , N, 
can be written as 

n eS (J 2 , J 4 , N; a^) = -J cff J2 ^ (m M m) > (3- 1 ) 

e£E* 

where the effective nearest-neighbor coupling reads 

Jeff =J 2 + J 4 J2 °l n) °j n) ■ ( 3 - 2 ) 

When all pairs of spins and for n ^ m are equal, the effective interaction 
J c ff = J 2 + {N — 1) J4 = J is finite. Otherwise, one has the coupling J cfT = J 2 + (N — 1 — 
2k) J 4 = J — 2k J 4 — ■> 00. Thus, the effective coupling is no longer translation-invariant. 
Nevertheless, the effective coupling J e g is always ferromagnetic. In this case, when the SW 

The word 'readily' here might seem misleading, since in some recent literature, the very inefficient 
Metropolis algorithm was still applied to simulate the standard Ashkin- Teller model in three dimensions, 
such as in Ref. [19]. 
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cluster algorithm is applied to update the spins, no frustration phenomena occur as in 
cluster simulation of the random Ising model. Therefore, it is plausible that the efficiency 
of such a cluster algorithm is not influenced too much by the inhomogeneous effective 
couplings. 

A Swendsen-Wang step for updating a given Ising variable can then be written 

as 

• Step 1: Place occupied bonds. Given a spin configuration a, for each edge 
e G E* of Q* , place an occupied bond with probability p 




1 — e 2J if <T- m ' 1 = <7 ? - m ' > and = for all n ^ m 

1 if oj"^ = of 1 ^ 1 and there is at least one pair ^ (3.3) 

n . r (m) I (m) 



• Step 2: Construct clusters. Two vertices that are connected through a chain of 
occupied bonds are said to be in the same cluster. On the basis of occupied bonds, 
the whole set of vertices V* of graph Q* is decomposed into a number of clusters C. 

• Step 3: Update spins a^ m \ For each cluster of C, randomly assign the value of 
all spins on its vertices to be +1 or — 1 with equal probability. 

This completes a Swendsen-Wang cluster step for spin variable a^ m K Since all other spin 
variables with n ^ m are kept fixed, this cluster method is called an embedding 
cluster algorithm. We shall refer to this method as the embedding SW algorithm for the 
IAT model. 

A valid embedding algorithm for the A^-color AT model must involve the update of 
all colors of Ising spin variables a. This can be done sequentially or randomly. 

4 Cluster simulation of the O(N) loop model 

Owing to the exact mapping between the loop and the IAT model, for integer > 
1, the embedding method described in Sec. II is already a valid cluster algorithm for 
the O(N) loop model on graph Q . In particular, for N = 1, this cluster method just 
reduces to the standard Swendsen-Wang cluster algorithm for the Ising model on graph 
Q*. The procedure in Sec. II was already applied to simulate the N = 2 IAT model on 
the honeycomb, square, and triangular lattices [21]. In comparison with the Metropolis 
method, the critical slowing down is considerably suppressed, but still observable. At the 
critical point of the triangular IAT model, which corresponds both to the critical Baxter- 
Wu model and to the 0(2) loop model on the honeycomb lattice, it was found [21] that 



8 



the dynamical exponent z is about 1.18. In this section, we shall show that the cluster 
algorithm in Sec. II can be further improved such that it can simulate the loop model 
with noninteger N > 1 and that critical slowing down hardly exists. 

Let us consider the Ising spin s = Ylni 17 ^- As mentioned earlier, any loop con- 
figuration of the O(N) loop model on Q is also a low-temperature graph of the s-spin 
configurations on Q*. But information about the colors of loops cannot be encoded in the 
s-spin configurations. 

In terms of loop configurations, the simulation of a^-spin variable in the embedding 
SW procedure in Sec. II is equivalently to update the loops of color m while keeping all 
other loops of color n ^ m unchanged. 

Our first task is then to rewrite the embedding SW procedure in Sec. II by using 
spin variable s and loop colors as the dynamical variables, instead of using o r( - m - ) for 
m — 1, 2, • • • , N . This leads to the following four steps: 

• Step 1: Place occupied bonds. Given a spin configuration s on Q* and the 

color information for loops on Q, for each edge e G E*, place an occupied bond with 
probability p 



• Step 2: Construct clusters on the basis of occupied bonds. Note that the s spins 
on the vertices of a given cluster can have different signs. 

• Step 3: Update spins s. For each cluster, flip all spins on its vertices with 
probability 1/2. 

• Step 4: Update loops. Consider the low-temperature graph of the s spin config- 
uration, let the colors of the existing loops remain unchanged, and assign the color 
of all newly generated loops to be m. 

Do Steps 1-4 for each of the N colors of loops. Step 4 is redundant for the procedure in 
Sec. II, since the colors of loops are already encoded in the spin configuration with 



For any satisfied loop configuration of the O(N) loop model, the color of each loop can 
randomly take any value of N colors with equal probability. Therefore, instead of keeping 
all loop colors n ^ m unchanged, one can randomly reassign the color of each loop in 





if Si 7^ Sj and e crosses a loop of color n ^ m 
if Si 7^ Sj and e crosses a loop of color m . 



(4.1) 



m = 1,2, ••• ,N. 
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Step 4, irrespective of the existing color information. On this basis, one can reformulate 
the embedding SW procedure as 



• Step 1: Construct loops and assign loop colors. Given a spin configuration s 
on £/*, construct the loop configuration on Q. For each loop, assign its color to be 1 
with probability 1/N and to be with probability 1 — 1/N. We say loops of color 
1 and to be 'active' and 'inactive', respectively. 

• Step 2: Update the s-spin configuration. This is done by the Swendsen- Wang- 
type procedure, in which the bond-occupation probability is 



A cluster simulation of the O(N) loop model can then be achieved by repeating the above 
two steps. Note that iV is no longer required to be an integer; it is sufficient to demand 
that N > 1, since 1/N should not be bigger than 1. 

In summary, in comparison with the embedding SW procedure described in Sec. II, 
the last version of the algorithm has two prominent features. First, N is no longer required 
to be an integer. Second, after a completion of the update of Ising spin variable, the loop 
colors are randomly reassigned, irrespective of the existing colors. This is possible because 
all the N colors of spin variables = 1, 2, • • • , N) are symmetric to each other. In 

other words, the symmetry among the iV-color spin variables is used in the last version 
of the embedding cluster algorithm. If this symmetry is broken in some way, e.g., one 
can let couplings J2 depend on spin colors, the last version is then not applicable. But a 
slightly modified version of the procedure in Sec. II still works. 

Because of usage of the symmetry among the iV-color spin variables, one would expect 
that modes of the critical slowing down, associated with this symmetry, are significantly 
suppressed in the last version of the embedding SW algorithm. Thus, one expects that it 
is more efficient than the procedure in Sec. II. This will be numerically confirmed later. 
Remark: For a plane graph Q of degrees d > 3, the loop assignment based on the 
s spin configuration is not unique, since loops are allowed to share a common vertex. 
Nevertheless, this should not affect the validity of the present embedding SW cluster 
algorithm. 

We conclude this section by pointing out that our embedding SW cluster algorithm 
for the O(N) loop model is of a similar spirit to the Chayes-Machta algorithm [22] for the 
g-state Potts model with real q > 1. 





if Si 7^ Sj and e crosses a loop of color 
if Si 7^ Sj and e crosses a loop of color 1 . 



(4.2) 
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5 The Honeycomb O(N) loop model 



The Honeycomb lattice is of degree 3, and thus the loops of the O(N) loop model never 
intersect. In this case, it has been known that the loop model can be exactly mapped on 
an O(N) spin model with partition sum [6-8] 

Z(w,N)= j\{dS k Hil + wSi-Sj), (5.1) 

k (ij)eE 

— * 

where E is the edge set of the honeycomb lattice. Here, S represents a N- dimensional unit 
vector - namely \S\ = 1 5 . This model has been under intensive investigation in the past 
two decades, and much exact information is available. For N < 2, the model undergoes 
a second-order phase transition; for N = 2, the transition is of infinite order-i.e., the 
so-called Kasteleyn-Thouless (KT) transition; for N > 2, it displays a lattice-gas-type 
transition. The renormalization exponents for critical O(N) loop systems for iV < 2 is a 
function of N. For N > 2, it is expected that the transition is in the same universality class 
of Baxter's hard-hexagon lattice gas model on the triangular lattice (an exact mapping 
between the N — > oo loop model and the hard- hexagon lattice gas can be established). 
The critical frontier for N < 2 is exactly located at [7, 8] 

w c (N) = (2 + y/2-Nj V2 . (5.2) 

It was observed [7,8] that a critical O(N) model corresponds with a tricritical q = N 2 - 
state Potts model. From exact calculations, conformal field theory [23], and Coulomb gas 
theory [8], exact values of renormalization exponents can be obtained for the tricritical 
Potts model as 



Vti = 3 

Vt2 = 4 

Vt3 = 5 

Vhi = 2 



6 

9 
16 

9 
30 

9 

(6-g)(g-2) 
89 



V» - 2- (1Q ~f + 2) , (5.3) 



5 In Ref. [7], the 0(N) spin model is defined such that |5| = VN and the factor w is replaced by X/N 
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where the Coulomb-gas couplings g is related to q as 



(7 = 4 + 



2 arccos 




(5.4) 



7T 



The symbols y t \, yt2, and y t 3 represent exponents of the leading, subleading, and next- 
subleading thermal scaling fields, respectively, and y^i and y^ are for the magnetic scaling 
fields. Therefore, exact values of critical exponents for the O(N) model described by 
Eq. dsnn function of N are also known. 

6 Simulation and Sampled quantities 

We applied the embedding SW cluster algorithm in Sec. Ill to simulate the 0(N) 
loop model on the honeycomb lattice; the Ising variable s is defined on the vertices of the 
dual lattice-i.e., the triangular lattice. The sampled quantities can be classified into three 
types; they are associated with loop configurations, Ising spin variables, and distributions 
of clusters C formed in the SW step. 

The first type includes the total lengths of loops and of the number of loops, as 
normalized by the volume of the system V and denoted by pi and p c , respectively. Here, 
the volume V is the total number of vertices on the corresponding triangular lattice. The 
second-moment fluctuations of pi are also sampled as C = V((p 2 ) — (pi) 2 ), where symbol 
( ) represents the statistical average. 

Taking into account the mapping between the loop model and the IAT model, we also 
measured both the total lengths and the number of 'active' loops, denoted by p\ a and 
p ca ; the fluctuations of pi a , denoted as C a , were also measured. The subscript 'a' means 
'active'. 

The second type of quantities concern the Ising variables associated with loop configu- 
rations. A natural quantity to sample is the magnetic susceptibility \ = V(m 2 ), where m 
is the magnetization density for spin variable s. After each embedding SW step, we also 
reconstructed Ising spins from these active loops. The associated susceptibility Xa 
was then sampled. In addition, we consider the size distribution of domains enclosed by 
those active loops. These domains can be obtained by placing occupied bonds between all 
equal nearest-neighbor Ising spins. They are named the Ising clusters in the standard 
Ising model. We measured the second moment of these Ising clusters 




(6.1) 
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where at is the number of spins in the ith domain. 

The last type of quantities contains information for the embedding SW simulations. 
They include the number of occupied bonds per spin n&, the number of clusters per spin 
n c , and the second and fourth moment of cluster sizes 

i 

where k = 2 or 4, and q is the size of the ith cluster formed in the embedding SW 
simulations. 

In the conventional SW simulation of the Ising model (or more generally of the Potts 
model), quantities W2 and W4 are exactly related to the second and the fourth moment of 
the magnetization density. The bond-number density is also exactly related to the energy 
density in the Ising model. 

On the basis of the above sampled quantities, we sampled ratios 

n-M! n {x{1))2 a n ^ 

^ ~ (x 2 ) ' ^ ~ (X {1) 2 ) ' ~ (3W| - 2W 4 > ' [ ' 

At criticality, these amplitudes ratios are dimensionless and universal. They are known 
to be very useful in locating critical points in Monte Carlo studies of statistical systems 
undergoing phase transitions. For the SW simulation of the Ising model, quantity Q w 
reduces to Q defined on the basis of magnetic susceptibility. 



6.1 Test of the algorithm 

For a test of the algorithm, we performed simulations for iV = 1.5 and N = 1.9. The 
system sizes took 12 values in range 8 < L < 256, where L is the linear system size of the 
triangular lattice. Periodic boundary conditions were applied. Measurements were taken 
after every Swendsen-Wang step. 

For N = 1.5, Eq. JOD yields the critical point at J c = -(lnio c )/2 = 0.24897- •■. 
Simulations were performed in range 0.245 < J < 0.253. About 6 x 10 6 samples were 
taken for each size L. Parts of the Q data are shown in Fig. El The intersections of the 
Q data for different system sizes rapidly converge to the expected value J c = 0.24897 
This implies the validity of the embedding cluster algorithm for the O(N) loop model 
with noninteger N. 
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Figure 2: Binder ratio Q for iV = 1.5. The data points +, x, □, 0> A, 0, and * represent 
system sizes L = 8,16,32,48,80,160, and 256, respectively. The error bars of the data 
are smaller than the point sizes. The lines, which simply connect data points for each L, 
are just for illustration. 

According the least-squares criterion, we fitted the Q data by 

m 

Q(J, L) = Q C + J^iJ - Jc) k L kyt + kL yi + hL^ + ■■■ , (6.4) 

k=l 

where m > 1 is an integer. The terms with amplitudes \ and b\ describe finite-size 
corrections. The term with exponent is supposed to arise from the leading irrelevant 
scaling field. Thus, the value of is given by y t3 in Eq. (|5.3|) . For N = 1.5, Eq. (J5.3|) 
yields = y t3 = —1.097296. Other finite-size corrections can arise from various sources, 
such as from the subleading irrelevant scaling field, from the analytical part of the free 
energy /, and from the second derivative of / with respect to the leading irrelevant field. 
It is not a prior clear how many or which such correction terms can be observed in the 
numerical data. Thus, one should make various fits by taking into account all possible 
sources of finite-size corrections. Final estimates of interested quantities are then obtained 
by comparing the results of various fits. We let exponent yt and yi to be fitted by the 
numerical data, and fixed y\ at —2 for simplicity. All the Q data in range 8 < L < 256 
can be successfully described by Eq. (J6.4J) . The results are y t = 0.749(4), yi = —1.4(3), 
J c = 0.24896(2), and Q c = 0.7429(8). The value of y t is consistent with the exact value 
Dt2 — 0.7481 ■ • • in Eq. (J5.3|) . and that of yi agrees with —1.097 • ■ ■ . The estimated critical 
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Figure 3: Binder ratio Q w for iV = 1.5. The data points +, x, □, 0> A> (), and * 
represent system sizes L — 8, 16, 32, 48, 80, 160, and 256, respectively, he error bars of the 
data are smaller than the point sizes. The lines, which simply connect data points for 
each L, are just for illustration purpose. 

point is also in excellent agreement with the exact prediction 0.24897- • ■ . 

Next, we plotted parts of the Q w data in Fig. EI The intersections of the Q w data 
reflect the percolation threshold of clusters formed in the embedding SW procedure. In 
cluster simulations of statistical systems, the efficiency can be reflected by the average size 
of the formed clusters. The efficiency will be limited if the average size is either too small 
or too big: little change of configurations is made for the former, and large amount of 
effort has to be carried out for the latter (probably it also accompanies by little change). 
Ideally, the percolation threshold of clusters formed in Monte Carlo simulations should 
coincide or be very close to the phase transition of statistical systems. This is indeed the 
case in the present embedding SW algorithm for the O(N) loop model, as reflected by 
Fig. El 

The fits of the Q W (L) data by Eq. (Q yields y t = 0.756(6), y { = -1.2(3), and 
J c = 0.24900(3). The estimated percolation threshold is consistent with the thermal 
critical point J c = 0.24897 

We also simulated the N = 1.9 loop model. Equation (|5.2|) predicts the critical point 
at J c = 0.20998- ■•. Simulations were performed in range 0.206 < J < 0.214, and 
system sizes took 12 values in range 8 < L < 256. Parts of the Q data are shown in 
Fig. El As expected, the intersections of the Q data for different sizes converge rapidly to 
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Figure 4: Binder ratio Q for N = 1.9. The data points +, x, □, 0> A, (), and * represent 
system sizes L — 8, 16, 32, 48, 80, 160, and 256, respectively, he error bars of the data are 
smaller than the point sizes. The lines, which simply connect data points for each L, are 
just for illustration purpose. 

J c = 0.20998 • ■ ■ . The fits of the Q data by Eq. (Q yields y t = 0.373(6), y { = -1.2(2), 
J c = 0.21000(8) « 0.20998- ••, and Q c = 0.586(1). The value of y t is more or less 
consistent with that of y t 2 = 0.367- ••. But that of yi does not agree with that of 
y t3 = — 1.81 ■ ■ • . In this case, we expect that the dominant finite-size corrections are not 
described by exponent y t 3- 

A plausible scenario for the dominant finite-size corrections might be the following. 
As mentioned earlier, the mapping between the O(N) model on graph Q and the I AT 
model on the dual graph Q* is only exact up to boundary effects. Even though these 
boundary effects are expected to vanish in the thermodynamic limit L — > oo, they cannot 
be neglected for finite L. According to a simple argument, these boundary effects vanish 
as a function of 1/L. Indeed, exponent —1 is consistent with the numerical values of yi 
for N = 1.5 and N = 1.9. 

7 Simulation at criticality 

For iV = 1.25, 1.5, 1.75, and 2, we performed extensive simulations right at the exactly 
predicted critical points given by Eq. (|5.2jl . Periodic boundary conditions were applied, 
and the system sizes took 18 values in range 4 < L < 1024. Samples were taken after 
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every embedding SW step. The number of samples is 10 for L < 256, and 2 x 10 6 for 
L > 256 (Wenan and Henk: Monte Carlo data are not complete yet). 

For later convenience, we list in Table. El the exact values of critical exponents for 
N = 1, 1.25, 1.50, 1.75, and 2, as predicted by Eq. ()5.3|) ; also included are the exact values 
of critical points given by Eq. ()5.2|) . 



N 


yn 


Vt2 


Vt3 


Vhi 


Vh2 


Jc 


1 


1.875 


1 


-0.625 


1.9479 • • 


1.1979- • • 


0.274653 • • • 


1.25 


1.8327 • • 


0.8873 •• 


• -0.8361 •• 


1.9343 • • 


1.1562 - • - 


0.263231- •• 


1.50 


1.7805 •■ 


0.7481 ■ ■ 


• -1.0972 •• 


1.9198 •• 


1.1069- •• 


0.248970 ■ ■ ■ 


1.75 


1.7078 •• 


0.5542 •• 


• -1.4607 •• 


1.9034 •• 


1.0420- •• 


0.229072 • ■ ■ 


2.00 


1.5 





-2.5 


1.875 


0.875 


0.173286- ■■ 



Table 2: Exact values of critical points and exponents for N = 1, 1.25, 1.75, 2, as predicted 
by Eqs. (Q and 



7.1 N = 2: KT point 

The critical N = 2 loop model on the honeycomb lattice is a very special case. It can 
be mapped both onto the standard N = 2 IAT model on the triangular lattice and the 
Baxter- Wu model with three-spin interactions on the triangular lattice. The latter has 
been exactly solved by various methods, and is believed to be in the same universality class 
as the 4-state Potts model. Further, it is known that the N = 2 loop model undergoes a 
KT transition, which is an infinite-order phase transition. 

Energy-like quantities. The energy-like quantities include the density of the length 
of loops pi and of the loop number p c . The exact information for the critical properties of 
the 0(2) model tells one that, at criticality, the finite-size behavior of pi and p c is governed 
by the sub leading thermal exponent yti in Eq. (|5.3|) . From the analysis of Monte Carlo 
data, as shown later, it turns out that the finite-size scaling behavior of bond-occupation 
density n&, a quantity associated with the embedding SW step, is also described by y t 2- 
Instead of being listed in tables, the pi and p c data are shown in Fig. El the vertical axis 
shows the data for pi(L) — pi , p c (L) — p c0 , and rib(L) — ribo, where constants p;o, p c o, and 
ribQ were obtained from the numerical fits. 

We made a least-squares fit of the Monte Carlo data for energy-density-like quantities 
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Figure 5: Quantities pi, p c , and n b for N = 2. In order to display the finite-size depen- 
dence for different quantities, the analytical-background contributions in these quantities 
have been subtracted. Namely, the data in the vertical axis are pi(L) — p i0 , p c (L) — p c0 , 
and n b (L) - n b0 , where p m = 1.114835(4), p c0 = 0.0574189(7), and n m = 1.10956(2) were 
obtained from fits by Eq. (J7.1j) . The exponent —2 of L in the horizontal axis is y t 2 — 2. 



to the general ansatz 



The symbols E and E should be replaced by specific quantities in the fits. For instance, 
depending on which quantity is going to be analyzed, E can be pi, p c , and n b , and E can 
be pi 0} p c0 , and n b0 . 

In Eq. (17. 1|) . the term with bo describes boundary effects arising from the mapping 
between the Ising-spin configuration and the loop representation. The term with bi arises 
from the least irrelevant scaling field. However, the value of yi is not very clear. According 
to conformal field theory, if an operator with exponent y exists, a sequence of exponents 
y — 1, y — 2, • • • , can in principle also exist. Therefore, yi can be y t 3 or y t 2 — 1 etc. In the 
analysis of our Monte Carlo data, we have tried to take into account all possible sources 
of finite-size corrections. 

The fits yield y t = 0.1(2) and p m = 1.114835(4) for the Pl (L) data, y t = 0.004(8) and 
p c0 = 0.057418(7) for the p c (L) data, and y t = -0.5(10) and n b0 = 1.10956(2) for the 
n b (L) data. The values of y t are all in good agreement with the prediction y t2 = 0. 

Within the sampled quantities, the density of the lengths and of the number of active 



E(L) = E + L yt - d 



(a + boL- 1 + b t L* + •••)• 



(7.1) 
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Figure 6: Specific-heat-like quantity C for N = 2. The exponent —2 of L in the horizontal 
axis is equal to 2y t2 — 2. 

loops, pi a and p ci , are also energy-like quantities. Since the color of each loop can randomly 
take one of N colors, the statistical mean value of pi a {L) is just p c (L)/N for any size L; 
the same applies to quantity p c . Therefore, we do not need to analyze them. 

Specific-heat-like quantities. Quantities C and C a account for the fluctuations of 
pi and pi a , respectively. The C and C a data are respectively shown in Figs. El and These 
suggest that the finite-size behavior of C is governed by exponent 2y t2 — 2 = —2, while 
that of C a is governed by 2y n — 2 = 1. 

An observation is that the behavior of the magnetic susceptibility x-, f° r t ne s = Y[ 
for m = 1, 2, • • • , N, is also described by exponent 2y t i —2 = 1. Actually, for the present 
case N = 2, because of the symmetry between two Ising variables and for the IAT 
model, it can be shown that, in the thermodynamic limit L — >• oo, C a (L) and x{L) are 
equivalent to each other (Wenan Henk: this statement has not been carefully checked). 

The data for specific-heat-like quantities were fitted by 

C(L) = C + L 2 ^-\a + boL- 1 + + ■••), (7-2) 

where symbols C and Cq should be replaced by specific quantities in the fits. 

The fits of the C(L) data yield y t = —0.01(3) = y t2 and Cq = 1.7867(8), and those of 
C a give y t = 1.5002(4) = y tl and C a0 = -0.46(3). 

Susceptibility-like quantities. Quantity \ for the s Ising variable has been dis- 
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Figure 7: Quantities C a and x f° r N = 2. The exponent 1 of L in the horizontal axis is 
equal to 2y t \ — 2. In this scale, the C a (L) and x{L) data collapse for each L. 

cussed earlier. The fits of the x{L) data also yield y t = 1.5002(4) = y n . Other 
susceptibility-like quantities include Xa for spin S 2a for the Ising clusters of spin 
variable a^ 1 ', and W 2 for the clusters formed in the embedding SW step. The data of 
these quantities are shown in Figs. |H1 and El 

The data for susceptibility-like quantities were fitted by 

K(L) = /C + L 2 ^-\a + boL- 1 + b t L^ + •••), (7.3) 

where symbols /C and /Co should be replaced by specific quantities in the fits. Exponent 
yh is left to be fitted by numerical data. 

The fits of the Xa(L) and W 2 {L) data by Eq. fl£3|) yield y h = 1.8751(1) and y h = 
1.8751(2), respectively. Both results are in good agreement with the exact value of yti = 
15/8 = 1.875, as given in Table El This implies that the percolation threshold of the 
clusters formed in the embedding SW step coincides with the thermal critical point J C (N = 
2) = In 2/4. 

The fit of S 2a yield yh = 1.9444(1). There is no prediction of the exact value for this 
exponent. 

Binder ratios. In the Monte Carlo simulations, we also samped dimensionless am- 
plitude ratios Q and Q a , defined by Eq. ()6.3jl . The Q and Q a data are shown in Fig. [TU1 
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Figure 8: Susceptibility-like quantities Xa/L 2 and W 2 for N = 2. The exponent —1/8 of 
L in the horizontal axis is equal to 2y hl — 4. 
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Figure 9: Quantity f° r 1 
obtained from the fits. 
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The data for amplitude ratios were fitted by 

Q(L) = Qo + + ■■■ , (7.4) 

where exponent yi is left to be determined by numerical data. The fits of the Q(L) data 
yield Q Q = 0.4774(4) and y, t = -1.03(6), and those of Q a give Q a0 = 0.6927(4) and 
yi = —0.99(3). Both estimates of yi agree with integer —1, which can be explained as 
Vt2 ~ 1. 

7.2 N = 1.25, 1.50, and 1.75 

The analysis of the Monte Carlo data for the honeycomb O(N) model with iV = 
1.25, 1.50, and 1.75 follows an analogous procedure as that for N = 2. Thus, the details 
of the fitting procedure are skipped here. The fitting results are listed in Tables El and |U 

For N = 2, quantities associated with those active loops or the spin variable a^ 1 ' have 
clear physical meaning, since they appear naturally in the N = 2 IAT model described 
by Eq. (|2.2j) . In this case, the scaling behavior of specific- heat-like quantity C a and sus- 
ceptibility Xa can be predicted from the exact results for the Baxter- Wu model or for the 
triangular IAT model. The numerical results fit well with the predictions. Nevertheless, 
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N 




Pc 


n b 


C 


c a 


1.25 


y 


0.887(2) 


0.884(5) 


0.885(5) 


0.885(4) 


1.378(3) 




y (exact) 


0.887- •• 


0.887- •• 


0.887- •• 


0.887- ■• 


- 




r 


0.61297(2) 


0.040207(2) 


1.09962(2) 


8.9(1) 


2.7(2) 


1.5 


y 


0.745(3) 


0.74(1) 


0.75(2) 


0.74(1) 


1.406(2) 




y (exact) 


0.748 ■ ■ • 


0.748 ■ • • 


0.748- ■• 


0.748- ■■ 


- 




r 


0.72952(2) 


0.046245(2) 


1.13371(2) 


4.8(2) 


1.21(4) 


1.75 


y 


0.54(1) 


0.54(1) 


0.50(5) 


0.57(4) 


1.432(3) 




y (exact) 


0.554-.. 


0.554.-. 


0.554--- 


0.554- •• 


- 




r 


0.86058(1) 


0.051884(2) 


1.15514(2) 


3.160(6) 


0.7(9) 


2.00 


y 


0.1(2) 


0.004(8) 


-0.5(10) 


-0.01(3) 


1.5002(4) 




y (exact) 














1.5 




r 


1.114835(4) 


0.057418(7) 


1.10956(2) 


1.7867(8) 


-0.46(3) 



Table 3: Fitting results for energy-associated quantities. Symbol r represents those 
analytical contributions, arising from the regular part of the free energy. Also included 
are the predicted values for exponent y to be fitted. Symbol "-" means that no prediction 
exists. 



no exact results seem to exist for the finite-size behavior of the geometric quantity S2 a , 
the second-moment of the Ising clusters for spin 

For noninteger N, however, the symmetry between the active and the inactive loops 
breaks. To our knowledge, there are thus far no prediction for exponents governing the 
scaling behavior of quantities C a and %ai as listed in Tables HU and |H On the other hand, 
the following scenario seems plausible. The estimated exponent y for C a in Table El is 
actually the mixed effect of y t \ and y t 2, while that for Xa is due to the effective mixture 
of yhi and yya- We did try to make fits for the C a and Xa by taking into account such a 
mixture effect, but the results seem worse than those by Eqs. (|7.2|) and (|7.3|) . 

An observation is that, independent of the A" value, the exponents for the scaling 
behavior of Xa and W2 are alway equivalent. This strongly suggests that the percolation 
threshold of clusters formed in the embedding SW step coincides with the thermal critical 
point for any N > 1. A rigorous proof is lacking. 
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N 


Y 

A 




Wo 




Q 


Qn 


1.25 


y 


1.8327(2) 


1.8806(3) 


1.8804(3) 


1.9491(2) 


-1.6(2) 


-0.56(6) 




y (exact) 


1.8327- •• 


- 


- 


- 


-1.09 - •• 


- 




Q 










0.8069(2) 


0.8245(7) 


1.5 


y 


1.7801(4) 


1.8843(2) 


1.8841(2) 


1.9499(1) 


-1.1(2) 


-0.60(4) 




y (exact) 


1.7805- •• 


- 


- 


- 


-1.09- •• 


- 




Q 










0.7426(4) 


0.7931(4) 


1.75 


y 


1.7079(4) 


1.8854(2) 


1.8854(1) 


1.9497(1) 


-1.6(2) 


-0.79(8) 




y (exact) 


1.7078- •• 


- 


- 


- 


-1.46- •• 


- 




Q 










0.6585(2) 


0.7628(4) 


2.00 


y 


1.5004(5) 


1.8751(1) 


1.8751(2) 


1.9444(1) 


-1.03(6) 


-0.99(3) 




y (exact) 


1.5 


1.875 


1.875 




-1? 


-1? 




Q 










0.4774(4) 


0.6927(4) 



Table 4: Fitting results for susceptibility-associated quantities. Symbol Q represents the 
universal values of those amplitude ratios at criticality. We also give the predicted values 
for exponent y to be fitted. The exponent y for amplitude ratios is for the dominant 
finite-size corrections. The question mark means that we are not very sure about the 
prediction of y^. 



8 Dynamical behavior 

During the Monte Carlo simulations, the value of every sampled quantity was saved 
on hard disk after every sweep. Statistical analyses were then performed on these data. 
Both static and dynamic information of sampled quantities can then be obtained. 

Let symbol / be a given quantity, the unnormalized autocorrelation function is then 
defined as 

Cf(t) = </(0)/(f)> - (f) 2 . (8.1) 
The normalized autocorrelation function is 

p f {t) = C f (t)/C f (0) . (8.2) 

Typically, p(t) decays exponentially (~ e~l*l/ T ) for large t, and the exponential autocor- 



24 



relation time is defined as 



W = lim sup r T-r, • ( 8 - 3 ) 
i-oo — In |p/(t)| 

In addition, we define the integrated autocorrelation time 

j oo ^ oo 

Tint,/ = 2 E Pf® = 2 + E^/W • ( 8 - 4 ) 

t=— oo t=l 

Here, the factor 1/2 is purely a matter of convention; it follows the definition in Ref. [24]. 
The integrated autocorrelation time controls the statistical error in Monte Carlo measure- 
ments of (/). More precisely, the sample mean has variance 

var(J) « i(2r m/J )CV(0) for n > r, (8.5) 
where n is the total length of Monte Carlo simulations. 

8.1 N = 2 

Parts of the numerical data for the normalized autocorrelation function are shown in 
Fig E] for susceptibility Xi i n Fig- El the density of loop lengths pi, and in Fig. IT31 for 
susceptibility Xa for spin variable The good collapse of the p data strongly suggests 
that the critical slowing down is absent in the embedding SW simulation of the 0(2) 
loop model on the honeycomb lattice. 

For quantity Xa, after a single Monte Carlo step, the autocorrelation function drops 
to a value smaller 0.1. This means that two subsequent samples are almost effectively 
independent. 

Except for very small t, the data lines in Figs llH IT2l andE3 are rather straight; the 
scattering phenomena at the right-hand side are due to statistical noise. Therefore, the 
autocorrelation function is almost a purely exponential function of time: p w e~ l l T . 

The integrated autocorrelation time Tint defined in Eq. (|8.4|) was also measured for all 
sampled observables. Apart from the constant factor 1/2, quantity Ti nt is just the area 
underline the p data line as a function of time t. Parts of the Ti nt data are shown in 
Fig.d 

We fitted the r int data for the sampled quantities by 

r int (L) = r + L z {b x + b 2 L- x + •••), (8.6) 
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Figure 11: Normalized autocorrelation function p for susceptibility x- The data points 
A, 0, and * are for system sizes L = 160, 256, 360, and 512, respectively. The scattering 
behavior at the right-hand side is due to statistical noise. For p x > 0.01, the p data for 
different system sizes collapse rather well in this scale. This suggests that the critical 
slowing down is absent in our embedding SW simulations of the 0(2) loop model on the 
honeycomb lattice. 
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Figure 12: Normalized autocorrelation function p for the density of loop lengths p\. The 
data points A, 0, and * are for system sizes L = 160, 256, 360, and 512, respectively. The 
scattering behavior at the right-hand side is due to statistical noise. For p x > 0.01, the 
p data for different system sizes collapse rather well in this scale. This suggests that the 
critical slowing down is absent in our embedding SW simulations of the 0(2) loop model 
on the honeycomb lattice. 
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Figure 13: Normalized autocorrelation function p for susceptibility \a corresponding Isin 
spin variable The data points A, (), and * are for system sizes L = 160,256,360, 
and 512, respectively The scattering behavior at the right-hand side is due to statistical 
noise. For p x > 0.01, the p data for different system sizes collapse rather well in this scale. 
This suggests that the critical slowing down is absent in our embedding SW simulations 
of the 0(2) loop model on the honeycomb lattice. 




Figure 14: Integrated autocorrelation time T int for N = 2. The data points □, 0> A, (), 
and * are for susceptibility x, susceptibility Xa for active spins, density of loop lengths 
pi, the second moment W 2 of the clusters formed in Monte Carlo steps, and the bond- 
occupation density rib, respectively. For the last two quantities, T int is very close to 0.5, 
which means that autocorrelation function pit) is essentially zero for t > 0. 
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Figure 15: Autocorrelation function p pi for N = 1.25. The data points A, 0, and * are 
for system sizes L = 160, 256, 360, and 512, respectively. The scattering behavior at the 
right-hand side is due to statistical noise. 

where the dynamic exponent z is left free to be determined by the numerical data. We 
assume that correction exponent —1, which appears in static quantities, also exists in the 
dynamical data. Nevertheless, it turns out that, for L > 8, all the r^t data for quantities 
in Fig. El can be fitted by Eq. ()8.6|) with 6 2 = 0. The results are shown in Table El 
These show that, for all quantities shown in Table the integrated autocorrelation time 
Ti n t converge rapidly to constants as L — > oo; the values of these constants are all rather 
small. Quantity \ has the largest autocorrelation time Tint, and pi has the next largest 

Tint- 

8.2 N = 1.25, 1.50 and 1.75 

For iV = 1.25, 1.50 and 1.75, the general behavior of p and r for quantities listed in 
Table El is similar as that for N = 2. Namely, the autocorrelation function p(t) is almost 
a purely exponential function of time t. For quantity Xa, after a single Monte Carlo step, 
the value of p drops to a significantly small value, which is about 0.27 for iV = 1.25. The 
quantities that have the largest two values of correlation time r int are \ and pf, pi has the 
largest value for iV = 1.25 and x f° r the others. Therefore, we only show the pi data for 
iV = 1.25 in Fig.lToland the x data for iV = 1.50 and 1.75 in Figs. ITBI and fTTl respectively. 

The Ti n t data for the quantities in Table El are also shown in Figs. HH1 EH1 and EH for 
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Figure 16: Autocorrelation function p x for iV = 1.50. The data points A, <0, and * are 
for system sizes L = 160, 256, 360, and 512, respectively. The scattering behavior at the 
right-hand side is due to statistical noise. 
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Figure 17: Autocorrelation function p x for N = 1.75. The data points A, 0, and * are 
for system sizes L = 160,256,360, and 512, respectively. The scattering behavior at the 
right-hand side is due to statistical noise. 
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Figure 18: Integrated autocorrelation time r int for N = 1.25. The data points □, 0> A, 
0, and * are for x, Xa, Pi, W 2 and n b , respectively. 

N = 1.25, 1.50, and 1.75, respectively. 

From Fig.^Jfor N = 1.25, it is not completely clear whether the critical slowing down 
is absent or not for quantities pi and x- Nevertheless, the approximately straight line of 
quantity pi suggests that, even the critical slowing down exists, the integrated correlation 
time t may diverge in a logarithmic format or in a power law of L with a very small 
exponent z. 

We fitted the r int data by Eq. (|8.6|) . For the T int)Pl data for iV = 1.25, it turns out 
they can still be described by Eq. (J8.6|) when a correction term with 6 2 is also included. 
For the other quantities and for the other values N, the r int data for L > 8 all can be 
well described by Eq. (18.6)1 with 6 = 0. The results are shown in Table El the value of b\ 
is 2.5(8) for T int>Pl with N = 1.25. 

Taking into account the approximate linearity of the T intyPl data for N = 1.25 in Fig. ll8[ 
we also fitted the r inttPl data by 

r int (L) =T + Ti\nL + b x L- x . (8.7) 

Indeed, all the data can be well described by Eq. (|8.7jl : we obtain r = 1.17(4), r$ = 
0.79(2), and bi = 0.8(2). 
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Figure 19: Integrated autocorrelation time T int for N = 1.25. The data points □, 0> A, 
0, and * are for x, Xa, Pi, W2 and rib, respectively. 




1000 



Figure 20: Integrated autocorrelation time r int for N = 1.25. The data points □, 0> A, 
0, and * are for x, Xa, Pi, W2 and n&, respectively. 
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Quantity 


N = 1.25 


N = 1.50 


N = 1.75 


N = 2.00 


Pi 


TO 


18(4) 


5.2(3) 


2.70(2) 


1.566(3) 




z 


-0.06(2) 


-0.17(2) 


-0.60(7) 


-1.7(1) 


X 




7.7(2) 


5.18(10) 


4.55(3) 


3.221(6) 




z 


-0.17(2) 


-0.43(3) 


-0.80(7) 


-2.3(3) 


Xa 


TO 


2.18(2) 


0.82(3) 


0.640(2) 


0.722(2) 




z 


-0.36(2) 


-0.28(4) 


-0.64(8) 


-0.94(6) 


n b 




0.499(4) 


0.498(4) 


0.500(1) 


0.5002(4) 




z 


-0.65(3) 


-0.6(1) 


-0.8(1) 


-1.2(2) 


w 2 




2.7(1) 


0.92(1) 


0.577(3) 


0.5036(8) 




z 


-0.23(3) 


-0.42(4) 


-0.54(6) 


-2.6(8) 



Table 5: Fitting results for the integrated autocorrelation time r^. 

9 Discussion 

By making use of the equivalence of the loop configurations of the O(N) loop model 
and the low-temperature graph of the Ashkin- Teller model in the infinite-coupling limit 
(IAT), we formulated an embedding Swendsen- Wang-type algorithm for the O(N) loop 
model with real value N > 1. For N — 1, this algorithm reduces to the conventional 
Swendsen- Wang method for the Ising model. With some modifications, an embedding 
Wolff-type method (single-cluster version) is readily available. 

We then applied our cluster algorithm to the O(N) loop model on the honeycomb 
lattice. The numerical data reveal the finite-size scaling behavior of several quantities. 
The associated exponents are confirmed to be those exact values predicted by the Coulomb 
gas theory and by conformal field theory. 

The dynamical data strongly imply that the embedding cluster algorithm suffers little 
from critical slowing down. This is somewhat impressive in the sense that the dynamic 
exponent z of the Swendsen- Wang type algorithm must satisfy the Li-Sokal bound [25]: 
z > a/u, unless it can be proved that the amplitude of terms with exponent a/u vanishes. 
The value of a/u = 2y tl —2 can be easily calculated from Eq. ()5.3|h which yields 1.6654 • • • , 
1.5610 ■ • • , 1.4156, and 1 for for A" = 1.25, 1.50, 1.75, and 2, respectively. The absence of 
a dynamic exponent z > av in the embedding SW simulations of the O(N) loop model 
implies that the amplitude for terms with exponent a/u indeed vanishes. We argue that 
this is because our embedding SW cluster algorithm has made use of the symmetry among 
the A" colors of Ising spins. 

Similar scenarios exist elsewhere. For instance, it can be proved that, in Metropolis 
simulations of the Ising model, the dynamic exponent z must satisfy z > j/u. However, 



32 



the dynamic exponent z of the Swendsen-Wang algorithm is much smaller than j/u. This 
is because the symmetry between the up- and down-pointing Ising spins is fully taken into 
account in the Swendsen-Wang algorithm. 

To have a better understanding of our above argument, let us consider another version 
of Swendsen- Wang-type cluster method for the IAT or the O(N) model with integer N, 
as described in Sec. II. This algorithm directly simulates each color of the spin variables 
a ( m ) j n the IAT model for N = 1, 2, • • • , N. In the language of loop configurations, the 
algorithm only updates loops in the same color; it does not interchange or reassign colors 
of loops. In other words, the symmetry among loops in different colors is not taken into 
account. Such a cluster algorithm has been applied to the critical N = 2 IAT model on the 
triangular lattice, which is equivalent to the critical 0(2) model on the honeycomb lattice 
or the Baxter- Wu model. It was found that the dynamic exponent z for the integrated 
correlation time is 1.18(2); indeed, z satisfies the Li-Sokal bound: z > ajv. 

Instead of the individual spin variables a^ m \ the final version of the embedding SW 
algorithm simulates the product variable s = Y[& ■ After each update of spin variable 
s, the colors of loops are reassigned randomly, irrespective of the existing colors of loops. 
In this sense, it is natural that dynamic exponent z > ajv vanishes. 

One would then expect that the subleading exponent a'/ v = 2y t2 — 2 serve as a lower 
bound for the dynamic exponent of our embedding SW cluster algorithm. 

The N = 1 loop model is just the Ising model. From universality, one has a'/ v = 
2y t2 — 2 = for the 0(1) model on any planar graph. On the honeycomb lattice, Eqs. (J5.3|) 
and (|5.4JI tell that a! /v is a monotonically decreasing function of q = N 2 . Therefore, one 
has a' jv < for N > 1. In this sense, our cluster algorithm still satisfies the Li-Sokal 
bound. Since the value of a'/u decreases as a function of N, it is also expected that, as 
N increases, the value of Tj nt decrease. This is consistent with r in Table El 

It is clear from Sec. Ill that the embedding cluster algorithm described in this work 
cannot be applied to the N < 1 case. Nevertheless, since all the loop configurations are 
the low-temperature graphs of the s spin variable, and vice versa. It seems that, for 
N <m 1, a reweighting modification of the Swendsen-Wang simulation of the Ising model 
can still be useful. Such a procedure can be described as 

• Step 1. For a spin configuration {s}, in which the number of loops is c, generate a 
new spin configuration {s'} by using the Swendsen-Wang algorithm. 

• Step 2. Derive the loop information for the new configuration {s'}, and calculate 
the loop number d . Accept the new configuration with probability N c '~ c , namely, 
set s <— s'. Otherwise, keep the old spin configuration s. 
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Repeating of these two steps forms a valid 'cluster' algorithm for N < 1. For the O(N) 
loop model on the honeycomb lattice, it turns out that, for iV > 0.8 and small system 
sizes L < 100, this algorithm works pretty well. 
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